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^^ Abstract 

^^ For a stationary reaction-diffusion system on a two dimensional domain we use the continuation 
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and bifurcation software pde2path to numerically calculate branches of fronts between different 



h^ patterns, and localized solution branches, for instance spots embedded in stripes and vice versa. 

^Z^ Some of these branches show a snaking behaviour in parameter space. We use the Ginzburg- 

^^ Landau reduction to approximate the locations of these branches by a Maxwell point for the 

VO associated Ginzburg-Landau system. 
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d 1 Introduction 

s 

I— I Homoclinic snaking refers to the back and forth oscillation in parameter space of a branch of 
stationary localized patterns for some pattern forming partial differential equation (PDE). Two 

^ standard models are the quadratic-cubic Swift-Hohenberg equation (SHe) 

m 

CN dtu = -(1 + A^)^ + Aix + lyu^ - ^x^ u = u{t, x) G M, x^VtdW^, (1) 



and the cubic-quintic SHe 

a^^x= -(1 + A^)^ + Aix + z/^x^ -^x^ i/ = i/(t, x) G M, xGf^cM'^, (2) 



with suitable boundary conditions if il 7^ R^, and where A G M is the linear instability parameter 
and z/ > 0. In both equations the trivial solution i/ = is stable for A < Ac := 0, where in the ID 
case d = 1 we have a pitchfork bifurcation of periodic solutions with period near 2ti. For z^ > i^o ^ 0, 
C^ with z/Q = \/27/38 for (M) and z/q = for fc), the bifurcation is subcritical and the periodic branch 

starts with unstable small solutions U- and turns around in a fold at A = Ao(z^) < to yield 0{1) 
amplitude stable periodic solutions. Thus, for Aq < A < Ac there is a bistable regime of the trivial 
solution and 0(1) amplitude rolls. 

In the simplest case the localized patterns then consist of ID rolls over a homogeneous back- 
ground ix = 0, and in each pair of turns in the snake the localized pattern grows by adding a pair 
of rolls on both sides, and this continues for ever over the infinite line. See, e.g., |BK06t IBKOTt 
lBKL+09] for seminal results in this setting, and |CK09^ IDMCKlT] for detailed analysis using a 
Ginzburg-Landau formalism and beyond all order asymptotics. ID snaking in finite domains has 
been studied in, e.g., |HK09^ iDawOQ ]. In 2D, by rotational invariance, depending on the domain 



and boundary conditions multiple patterns may bifurcate from i/ = at A = 0, for instance so 
called straight rolls, and hexagonal and rectangular spots. There are a number of studies of lo- 
calized patterns over two dimensional domains |LSAC08| lALB^lOj , often combining analysis and 
numerics, and in more complicated systems like fluid convection [LJBKTT], but all these essentially 
consider patterns embedded in a homogeneous background. Already in |Pom86j it is pointed out 
that "pinned" fronts connecting stripes and hexagonal spots may exist in reaction-diffusion systems 
as a codimension phenomenon in parameter space, i.e., for a whole interval of parameters, and 
similar ideas were also put forward non-rigorously in |MNT90j with the 2D quadratic-cubic SHe 
([2]) as an example, see in particular |MNT90[ Appendix C]. Such a pinned front is observed in 
|HMBD 95] for ^ by time integration, but so far no studies of branches of stationary solutions 
involving different 2D patterns seem available. 

Here we consider a standard model problem for predator (u) prey (v) reaction diffusion systems, 
namely 

dtU = DAU + N{U, A), N{U, A) = (7_^X'') ' ^^^ 

with U = {u^v){t,x,y) G M^, diffusion matrix D = (oS)' ^ fixed to d = 60, and bifurcation 
parameter A G M+. The reaction term A^ of ([3]) is a special version of the Schnakenberg |Sch79| 
and Selkov model |Sel68| . In particular we consider the stationary system 

DAU + N{U, A) = 0. (4) 

The unique spatially homogeneous solution of Q is it;* = (A, 1/A). We write Q in the form 

dtw = L(A)w + G{w), (5) 

with w = U-w'' and L(A, A) = J(A)+i?A, where J is the Jacobian of A^ in it;*. For k = (m, n) G M^ 
we have L(A, A)e^(^^+^^) = L(k, A)e^(^^+^2/)^ where L(k, A) = J(A) - Dk^, k := Vm^ + n^, and 

thus we also write L{k^X). The eigenvalues of L are given by /x±(k. A) = l~^±{k^X) — ^ ^ =^ 
^ ^'^^^^^) \ -detL(A:,A). Following jMur89[ Chapter 14] we find that in Ac = \/d\/3 - \/8 ^ 



3.2085 we have /i+(k, Ac) = for all vectors k G M^ of length kc = \/\/2 — 1 ^ 0.6436, and all other 
/i±(k, Ac) < 0, i.e., there is a Turing bifurcation at Ac, with critical wave vectors k G M^ in the 
circle |k| = k^ The most prominent Turing patterns near bifurcation are stripes and (hexagonal) 
spots and can be expanded as 

[/ = !(;* + 2 ( Acos(A:cx) + Scos ( y (-x^V^y\\ -{- B cos (^ (-x- VSy)) j $ + h.o.t. 

= !(;* + 2 ( Acos(A:cx) + 25cos ( — ^) cos (—VSyj j $ + h.o.t., (6) 

where $ G M^ is the critical eigenvector of L{kc, Ac), and h.o.t. stands for higher order terms. The 
amplitudes 2^4, 2B G M (where the factor 2 has been introduced for consistency with ©1) of the 
corresponding Turing pattern depend on A, with A = S = at bifurcation. 

In ^we numerically calculate Turing patterns for Q, including so called so called mixed modes, 
and, moreover, find some branches of stationary solutions which involve different patterns, namely 

(a) fronts between stripes and spots, and associated localized patterns, e.g., stripes localized in 
one direction on a background of hexagons, and vice versa; 



(b) fully 2D localized patches of hexagons over a homogeneous background. 

The patterns in (a) come in two regimes: one "hot" which is rather far from the primary bifurcation 
at Ac; one "cold", with A closer to Ac- Of these, only the "hot" branches show snaking behaviour in 
our numerical simulations, and we conjecture that the non-snaking of the "cold" branches (and also 
of the fully localized branches (b)) is due to the fact that subcriticality of the primary bifurcation is 
rather weak. In ^we relate our numerical results to arguments derived from the Ginzburg-Landau 
reduction. The main result is the calculation of three Ginzburg-Landau Maxwell points which give 
predictions for the A-ranges of the fronts and localized patterns in the full system Q. In ^we 
close with a brief discussion. 

The results are not specific to the model system Q but can be expected in any reaction diffusion 
system (over sufficiently large domains) with a bistability between different patterns which allow 
homo-or heteroclinic connections in the associated Ginzburg-Landau system. 

2 Numerical Results 
2.1 Stripes, spots, and beans 



We use the bifurcation and continuation software pde2path |UWR12| to numerically calculate 
patterns for Q, and their stability, over the domain ^ = [—Ix^lx] x [~lyi^y]i Ix — '^h^/kc^ ly — 
26127T / {\/3kc) , with Neumann boundary conditions, which was chosen to accomodate the basic 
stripe and spot patterns. Here h^h ^ N, and the slight "detuning" 5 ^ 1 is used to unfold the 
multiple bifurcation at A = Ac since pde2path currently only deals with simple bifurcations. 
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Figure 1: Bifurcation diagram and example solutions of Q obtained with pde2path over the "2 x 2" -domain 
^ = [—1x1 J^x] X [—J^yi J^y] with Neumann boundary conditions, lx=4:7r/kc^ ly=47r/{V3kc)^ kc = v a/2— 1^^0.6436. 
The branch horn are the homogeneous solutions, hs the hot (up) stripes, cs the cold (down) stripes, hh the 
hot spots, and ch the cold spots. For instance hsp30 stands for the 30*^ point of hs. hb and cb are mixed 
modes, which we call hot and cold beans, respectively. For hs, hh, and hb we plot the maximum of u, and 
the minimum for cs, ch, and cb. Stable and unstable parts of branches are represented by thick and thin 
lines, respectively. In particular, except for hhplO, hhp30 all displayed patterns are unstable. Only a small 
selection of bifurcation points is indicated by o. 



In Fig. [T] we use 6 = 0.99 and li — I2 — 2, which we cah a 2 x 2 domain as in both directions 
2 spots "fit". The main panel shows a part of the very rich bifurcation diagram (BD). Following 
biological terminology we classify the stripes (A 7^ 0, S = 0) into hot (also called "up") stripes (hs, 
A > 0) and cold (also called "down") stripes (cs, A < 0), which exhibit a maximum respectively a 
minimum at x = 0. They are 27r/A:c-periodic in the horizontal direction, bifurcate in a supercritical 
pitchfork, and are stable from A55 ^ 3.15 to Xge ^ 2.51, where here and henceforth the first index of 
the bifurcation parameter A stands for the branch and second for ending or beginning of stability. 
The (hexagonal) spots (A = S 7^ 0) can be classified into hot spots (hs, A^ B > 0) and cold spots 
(cs, A = S < 0), which have a maximum resp. a minimum in the center of every spot. They are 
47r/A:c-periodic in the horizontal and 47r/(\/3A:c)-periodic in the vertical direction. The spot branch 
bifurcates transcritically from the homogeneous branch at A = Ac, i.e., it starts in a saddle node or 
fold bifurcation Xch ^ 3.22, and the stability region of the cold spots begins at the fold and ends in 
Ace ^ 3.03. The hot spots are stable from A/^^ ^ 2.73 to A/^g ^ 0.98. 

Thus, there is a bistable range of cold stripes and cold spots for A G (Ace, ^sh)^ of hot stripes and 
hot spots for A G (A^e, A/^^), and of cold spots and the homogeneous solution for A G (Ac, Ac^). A 
branch of "skewed hexagon" or "mixed mode" solutions of the form ([g]) with A < B <{) bifurcates 
subcritically from the cold spot branch ch in Ace and terminates on the cold stripe branch cs in A55. 
We call this type of solutions cold heans (cb), cf. | Yan04j . There is also a branch hb of hot beans 
with A> B > {) which bifurcates subcritically from the hot stripe branch hs in A^e and terminates 
on the hot spot branch hh in A/^^. We always plot the patterns of u. If we have a hot pattern for 
u^ then i; is a cold pattern and vice versa. This follows from the predator-prey structure of the 
reaction term A^ of Q. 

Remark 2.1. For all these patterns the discrete symmetry S = S2T^/kc gi^^^i by the shift by 27r/A:c 
in the horizontal direction yields a solution as well. The stripes are of course invariant under 5, but 
for the hexagonal spots the shifts generate new branches 5hh and 5ch. These make the plotting 
of the BD a bit complicated graphically. Therefore, in Fig. [2] we repeat the BD from Fig. [l] in a 
schematic way, with ordinate ix(0, 0), i.e., u in the center of the domain, and add the branches of 
shifted hexagons. The bean branches hb and cb then both take part in a "loop" involving shifted 
beans and so called rectangles which again can be expanded as in ^ but with \A\ < \B\. This can 
be worked out on the level of amplitude equations, see ^ and it is also recovered by our numerics 
for the full system. The behaviour described below like snaking branches bifurcating from hot 
beans obviously transfers to branches related by 5, while on the rectangle branches some different 
behaviour occurs, which here we do not study, but see [WF13J . In the following we mainly focus 
on the hb and cb branches, and also mostly restrict to following just one direction at bifurcation. J 

Remark 2.2. The main regime of interest to us is A G [2.5, 3.22], and below A ^ 2.5 all branches 
plotted except hh are unstable. We plotted a somewhat larger bifurcation picture since, e.g., on 
the ch branch some interesting patterns occur. J 

Remark 2.3. pde2path uses the Matlab FEM pdetoolbox to discretize elliptic PDEs like 
including some error estimators and adaptive mesh-refinement, see |UWR12] for details. In Fig. 
we used a regular "base mesh" of 20.000 triangles which, e.g., on the beans branch is refined to 
about 60.000 triangles on average. Moreover, for all solutions calculated the mesh can be further 
refined to yield arbitrary small error estimates without visible changes to the solutions. Calculation 
time for the full BD in Fig. [T]on a quad core desktop PC is about 40 Minutes. 

The rather large numbers of triangles was needed mainly in order to avoid uncontrolled branch 
switching, and it is the large number of different branches which makes the continuation and 
bifurcation numerics of Q demanding. For instance, on the 2x2 domain the first 10 bifurcations 
from the homogeneous branch occur between A = Ac ^ 3.2085 and Aio = 3.1651. On larger 
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Figure 2: (a) Schematic BD with 2/(0,0) as ordinate. Additional branches compared to Fig. [T] are dashed, 
and 5* denotes the respective phase shifted pattern, i.e., Shh are the phase shifted hot hexagons, hb and cb 
both form part of a loop hs ^ hh ^ Sch -^ cs -^ ch ^ Shh -^ hs, where rl and r2 are called rectangles. 
(b)-(d) rl example solutions for which, in the expansion (|6|, 5 > and A changing from A>0 in (b) via 
A=0 in (c) to A<0 in (d). (e)-(f) r2 with B < and A changing from A<0 in (e) via A=0 in (f) to A>0 
in (g). All domains as in Fig. IT] 

domains, these bifurcation points collapse to Ac, e.g., over the 6x2 domain used in Fig. [4] below 
we have Aio ^ 3.2042. Similarly, on all branches shown in Fig. [l] there are very many bifurcation 
points, and thus we only plot a small part of the full BD. See jUWR12| on details how pde2path 
tries to avoid undesired branch-switching. The pde2path script to generate Fig. [l] is also included 
in the software as demo schnakenberg. J 



2.2 Fronts, localized patterns, pinning and snaking 

Counting from the hs branch, pde2path yields four bifurcation points hbbpl, hbbp2, hbbp3 and 
hbbp4 on the hot bean branch. One bifurcating branch connects hbbpl and hbbp4, and another 
branch connects hbbp2 and hbbp3. 

By doubling the horizontal length to a 4 x 2 domain, i.e., setting Ix = Sn/kc, we find eight 
bifurcation points hbbpl,..., hbbp8 on hb, see Fig. [3] The branches bifurcating in hbbpl and hbbp2 
are called hf (hot front) and hi (hot localized), respectively. Some example solutions on hf and hi 
are also presented, and the so called Ginzburg-Landau Maxell point A^. The hf branch connects 
hbbpl to hbbp8, and contains stationary fronts t/het from hot spots to hot stripes, while hi con- 
nects hbbp2 to hbbp7 and contains homoclinic solutions t/hom in the form of localized hot stripes 
embedded in hot spots. Both branches show a snaking behaviour in the BD, and we first discuss 
the hi branch in (a) which indicates the start of so called homoclinic snaking, which becomes more 
prominent if we further increase the domain size. 

The idea is that on hi, looking from hbbpT to hbbp2, during each cycle consisting of two folds, 
a further roll is added on both sides of the roll pattern localized in the middle over a background 
of hexagons. After every second saddle-node the branch contains stable solutions. In contrast to 
theory and also to ID problems over really large domains, the branch does not snake around a 
vertical line but in a slanted manner. This, and the fact that the branch does not directly bifurcate 
from the stripes but from the beans, are finite size effects, cf. [ BCROSi IBBKMOSi IHKOQi lDaw09] . 
During, e.g., the initial traverse from hbbpT to hi 172 the "near hexagon bean pattern" at the 
bifurcation point is reshaped to spots at the sides and stripes in the middle. The analogous 
reshaping to "near stripe beans" takes place between hl96 and hbbp2. By symmetry there is also 
a branch of localized hot spots embedded in hot stripe patterns (not shown). 



Remark 2.4. For, e.g., the ID cubic quintic SHe Q, an analogon of the bean branch is the branch 
r_ of smah amphtude unstable rolls that bifurcates subcritically from (i/. A) = (0, 0) and "connects" 
to the fold where the rolls become stable. In [Daw09] the appearance of bifurcation points on r_ is 
related to modulational instabihty of the amplitudes ^-(A) of the unstable rolls in the associated 
Ginzburg-Landau equation, derived by using a scaling of the "subcriticallity parameter" z/ in ([2]). 
In particular, this gives a lower bound on the domain size necessary for these bifurcations, which 
is linear in 1/z/, i.e., inversely proportional to the subcriticality parameter. Moreover, it explains 
why the bifurcations occur in connected pairs like (in our case) hbbpl - hbbp8, hbbp2 - hbbp7, and 
so on. 

At least the relation between subcriticality of the equation and necessary domain size for sec- 



ondary bifurcations also holds for our system, see Remark 2.5 We expect that calculations similar 
to those in |Daw09| can also be done in our case using the system of Ginzburg-Landau equations 
derived in ^ but naturally they will be more complicated. J 
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Figure 3: Bifurcation diagrams and example plots of fronts and localized patterns of Q obtained with 

/ \ 1/8 

pde2path over a 4 x 2 domain. Here ||ii|| = Hi^H^s *= ( t^ /^ |'^(^7 ^)|^d(x, ?/) J , which was chosen purely 
for graphical reasons, and A^ ~ 2.62 is the so called (Ginzburg-Landau) Maxwell point, see ^ (a) shows 
the branch of localized patterns hi connecting the second and the seventh bifurcation points of the hot 
bean branch. Bifurcation points are indicated by o, stable and unstable parts of branches by thick and thin 
lines, respectively. hbbp2 stands for the first bifurcation point on hb. (b) shows the branch hf of hot fronts 
connecting hbbpl and hbbp8. 



In contrast to the homoclinic snaking in (a) which at least for simpler models on infinite cylinders 



can be also explained and predicted theoretically, see ^ finite size and numerical effects should 
be regarded as essential for the "heteroclinic snaking" in (b). In fact, at least for the ID Swift- 
Hohenberg equation, beyond all order asymptotics [CK09J predict that for each A near the Maxwell 
point there are at most two pinned front solutions, again see ^ Thus we believe that the fronts in 
(b) should be seen as homoclinics by even extension of the solutions over the left boundary (spots 
on a background of rolls) or right boundary (stripes on a background of spots). 



(b) cf 15, cf40 and cl24, cl50 
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Figure 4: Bifurcation from the cold beans on a 12 x 2 domain, (a) cf branch connecting cbbpl with cbbp8 
in a monotonous way (no snaking), cold localized branch cl, and the cold Maxwell point (see Q. (b) some 
example solutions. 

On the cold bean branch cb there are no bifurcation points over a2x2or4x2 domain (see 
Remark 2.5 below for further discussion). However, Figj4[a) shows 4 of the 8 bifurcation points 
obtained on the cold bean branch over a 12 x 2 domain, namely cbbpl, cbbp2, cbbp7, cbbp8, 
the branch cf (cold front) connecting cbbpl and cbbp8, and the branch cl connecting cbbp2 and 
cbbp7, and consisting of spots on a stripe background. Panel (b) shows 4 example solutions. In, 
e.g., cl50 it can be seen that the localized spots are shifted by 2t{ jkc in x, such that here we really 



consider 5cb which connects cs to 5ch, cf. Remark |2.1| but for simplicity we omit the S. 

Similar to Fig.|3j after bifurcation from, e.g., cbbpl, in the traverse to cf 15 the solutions reshape 
into stripes on the left and spots on the right. The growth of the hexagon part then happens in a 
narrow A regime around A^ between cf 15 and cf70, while between cf70 and cbbpB the solutions 
are reshaped to spot-like beans. Similar remarks hold for the cl branch. 

Remark 2.5. An essential difference to Fig. [3] is that there is no snaking behaviour of the cf 
and cl branches. Following Remark 2.4 we believe that the need for a large domain to obtain 
bifurcations from the cb branch can be explained from the "weak subcriticality" of the system in 
this range. This means that the cold bistable range is rather narrow, compared to the hot bistable 
range. In other words, the difference in (Ginzburg-Landau-)energy between the ch amplitudes and 
the cs amplitudes is much smaller than between the hh and hs solutions in their bistable range, 
and this results in flatter fronts for the associated Ginzburg-Landau system. 

However, here increasing the domain gives bifurcation points on cb (e.g., 4 on a 6 x 2 domain), 
and branches connecting the first and the last, the second and the second to last, and so on, as in 
|Daw09| . but no snaking. Presently we do not know if this non-snaking is still due to numerical 
and finite size effect, but rather conjecture that it is not: Increasing the domain size, the steep part 
near A^ becomes steeper, but our numerics suggest that this process converges to a monotonous 



in A branch, //there is snaking in these branches over finite domains, then it is in a very narrow 
A interval, and one needs very large domains and fine discretizations. 

Similar effects are also observed in [HO05J for the ID cubic quintic SHe ^ without a scaling 
assumption on z/. For 0<z/ = (9(l) there are branches which connect points on r_ (unstable small 
amplitude rolls) near zero with points on r_ near the fold, but only for sufficiently large (fixed) v 
they become steep in parameter A and eventually start to snake in the bifurcation diagram {v > 1.5, 
which agrees with the numerics in [DMCKIH Figure 1] for ^ with v — l.Q and a rather narrow 
snake) . J 
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Figure 5: Bistable range between homogeneous solution and cold hexagons, (a) partial BD and the homo- 
geneous Maxwell point (see Q. Additional to quasi ID solutions we find fully 2D localized hexagon patches. 
The triangle in (d) indicates the computational domain with a = (327r/A:c,0) and b = (327r//cc, 327r/(>/3A:c), 
and about 90.000 triangles. We use Neumann BC on all three sides and for plotting first make an even 
extension over y = and then five rotations by 60 degrees, (e) quasi ID front between w — (A, 1/A) and 
spots on a 13 X 1 domain, A ~ 3.219. 



The third bistable we discuss is between the homogeneous solution it;* = (A, 1/A) and the cold 
spots for 3.21 ^ Ac < A < X^h ^ 3.22. Here, over sufficiently large domains, we have bifurcation 
points on the small amplitude (unstable) cold spot branch, which now corresponds to the r_ 



branches for ^ mentioned in Remarks |2.4| and |2.5[ and as above we can find quasi ID fronts 
between it;* and cold spots, see Fig. [7|e), and the associated ID locahzed patterns. However, 
as in the qcSHe ([l]) and the cqSHe Q, see [ ALB^lOj and the references therein, we can also 
calculate fully locahzed patches of hexagons, see Fig. [7|(a)-(d). In contrast to |ALB+10| . none of 
these branches shows snaking (in our numerics, over any of the domains we considered), and again 
we conjecture that this is due to the very weak subcriticality, or, in other words, the very weak 
difference in Ginzburg-Landau energy between it;* and the cold spots. 

3 Ginzburg-Landau reduction 

We now approximate the spots, stripes and mixed modes by the Landau formalism, and the various 
fronts and localized patterns by the Ginzburg-Landau formalism. In particular, using an energy 
argument and the so called Maxwell point for the Ginzburg-Landau system we find an approximate 
prediction where to find the hot snaking and the cold fronts in the bifurcation diagram for Q. 
Similar ideas are worked out much more deeply for the ID quadratic-cubic SHe ^ in [CK09J and 
the cubic-quintic SHe ^ in [DMCK ll], where, by augmenting the Ginzburg-Landau ansatz with 
beyond all order asymptotics, accurate bifurcation diagrams for homoclinic snaking were rigorously 
derived. 

In our case we need a system of Ginzburg-Landau equations, and our analysis is more formal 
since even a consistent derivation of the Ginburg-Landau system is difficult as we refrain from 
scaling assumptions for quadratic interactions but rather work with the numerical coefficients in 
Q. However, from the Ginzburg-Landau system we calculate the Maxwell point as a necessary 
condition for Ginzburg-Landau fronts, and the pinning argument from [Pom86] then suggest the 
existence of stationary fronts for Q. This approximation turns out to be qualitatively and at least 
in the cold regime also quantitatively correct, and thus it gives a lowest order approximation for 
the numerical solutions, although it cannot explain the snaking. 

3.1 Landau description of spots, stripes, and mixed modes 

To formally reduce Q to a Ginzburg-Landau system the idea is to treat x as an unbounded variable, 
while y G [~~Wh, ~^Wh] with Neumann boundary conditions, as in the numerics. At least close to 
Ac the most unstable modes are then ei$, 62$, and 63$, where ei = e^^^^, 62 = Qikc{-x+vsy)/2^ 
63 = Qikc{-x-VSy)/2^ ^^^ ^ ^ ^f^^^j ^ ]^2 jg ^Yie eigenvector of L{kc,\) to the eigenvalue /i+(A:c, A). 
First we consider slowly varying complex amplitudes Aj — Aj{t) of these modes, j — 1,2,3, i.e., 
our ansatz reads 

3 3 3 

w = ^AiCi^ + 2 ^ l^iP^o + 5]] A^ef(/)i + ^ Ai'Ajeie-jC^2 + cc. + h.o.t., (7) 

1=1 i=l i=l l<^<j<3 

where Aj means the complex conjugate of Aj and cc. stand for the complex conjugate of all 
preceding terms. The vectors 0o, 0i, and 02 are introduced to remove quadratic terms at wave 
vectors k with A: = 0, /c = 2kc and k = V^kc from the residual L{A)w + G{w). 
The calculations are best organized by writing ^ in the form 

dtw = L{A)w + B{w, w) + C{w, w, w), (8) 

where B and C are symmetric bilinear and trilinear forms. To remove terms of order AiAj from 
Q we need 

0o(A) = -2L(0, A)-iB($,¥), 0i(A) = -L{2k„ \)-'Bi^, $), 

02(A) = -2L(^/3fce,A)-lB($,¥). 



These terms arise due to quadratic interactions of the forms, e.g., B{Aiei^, Aiei^) = |Aip5($, $), 
B(Aiei$, Aiei$) = A2e2i^-^S($,$), and S(Aiei$, A^e^) = Ai]4^e^^(^^-^^)S($, $). Allthough 
$ G M^ in our case we keep the notation $ as this makes it easier to see where the respective terms 
come from. The matrices L(0,A), L(2fcc,A) and L(\/3A:c,A) are invertible at least for A not too 
far from Ac- From the Fredholm alternative we obtain the Landau ODE system as the solvability 
conditions for removing terms up to cubic order, namely 



at ei : dtAi = /i(Ai, A2, A3) 
at 62 : dtA2 = f2{Ai,A2,As) 
at 63 : dtAs^ f3{Ai,A2,A3) 



ciAs + C2A1A2 + cs\As\^As + 64^3(1^1^ + IA2I 



ciAi + C2A2A3 + C3IA1IM1 + C4Ai{\A2r ■ 
C1A2 + C2A1A3 + C3IA2PA2 + 64^2(1^1^ ■ 



(10) 



with ci(A) = M+(fcc, 0, A], C2(A) = 2(S($, $),$*), C3(A) = (3C($, $, $)+2S($, 0i)+25($, ^o), $*), 
and C4(A) = (6C($,$,$) + 2S($,(/)2) + 2S($,(/)o), $*). Here $*(A) is the adjoint eigenvector of 
L(kc,X) to the eigenvalue ii^{kc,X), normalized such that ($,$*) = 1. At ej, j = 1,2,3, we ob- 

J K84llPom86j . pl93l §IV,A,l,a(iii)], |DSSS03[ 



tain the complex conjugate equations. See, e.g. 

§2], |Hoy06j §5],^NG06nPis06j, which also explain that (10) is the generic form of the amplitude 



equations for hexagonal symmetry. The crucial step is the actual calculation of the coefficients. 

Remark 3.1. In Q, $ = ^(/cc. A) varies with A and is not fixed at A = Ac, which would be the 
more classical ansatz. Similarly, terms like Q and the coefficients C2, . . . , C4 are evaluated at A, not 
as usual at A = Ac- The formalisms are in fact equivalent for A ^ Ac with the differences hidden 
in the h.o.t. in Q. The reason why we always evaluate at A is that we want to use the formalism 
also for Ac — A = 0{1) and then expect (and find) better approximations with $ = $(A). On the 
other hand, for the comparison with the numerics we want to keep the wave vectors fixed and thus 
evaluate at fee, which approximately stays the most unstable wave number also for Ac— A=(9(l). J 

From the coefficients Cj(A), j = 1, . . . , 4, see Fig. [oFa), it follows that 



T± 



±J- 



C3 



and P± = — 



C2 



Cl 



2(C3 + 2C4) V 4(C3 + 2C4)2 C3 + 2C4 ' 



(11) 



are real in the Turing-unstable range A < Ac, respectively for A < AcLfoid ^ ^cb ^ 3.22. The triples 
(T+, 0, 0), (T_, 0, 0), (P+, P+, P+), and {P-,P-,P-) solve ^ and via generate hot stripes, cold 
stripes, hot spots, and cold spots, respectively. Mixed modes are obtained from setting A3 — A2 
and solving (10) for (^1,^42). The formulas, together with the signs of Cj(A), also explain the 



supercritical pitchfork bifurcation of stripes, and the transcritical bifurcation of the spots. Figures 
|6][b),(c) show the complete BD for stationary solutions of ( 10) in the invariant subspace Ai =: A G M 
and A2 — A3 —\ B ^ R, while (d),(e) compares the BD of stripes, spots and beans for (10) with the 
numerical BD from Fig. [l]for the full system ([1]). Though qualitatively correct down to A = 2.4, 
the approximation errors grow with | A — Ad, as expected. Additionally to the L^ like error already 
shown, for i G {stripes, spots} we define the relative L^ errors 



E 



\^mim,i\X^y) — Ui[X^\ 



E 



\u„ 



n^5 ^ 



(12) 



between the numerical solutions t/num and approximate solutions U = w"^ + w from ([7|, where Q^ 
is the discretization used to calculate the numerical solution from Fig. fl] It turns out that ^stripes 
stays rather small also for Ac — A = 0(1), while Zhex behaves worse, see Fig. [OJ^f), and also Remark 
14.11 for further comments. 
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Figure 6: (a) Landau coefficients, and Ginzburg-Landau coefficient cq, see ^3.2| (b),(c) Landau BD in the 
invariant subspace Ai =: A G M and A2 = As =: B G R, stability/instability indicated by thick/thin lines. 
(d),(e) Comparisons of numerical BD from Fig. [l] (indicated by * for stable points and + for unstable points) 
with Landau BD near the cold and hot bistable regimes, (f) relative L^ errors between hot numerical and 



Landau solutions defined in (12). 



Another notable discrepancy between the Landau BD and the (numerical) BD for Q is that in 
the former the (hot) stripes are stable all the way to the bifurcation point, while in the PDE this 
depends on the domain size. On a 1 x 1 domain we obtain exactly the same stability as displayed 
in Fig. m[b),(c), and the instability near bifurcation of the stripes already on a 2 x 2 domain as 
in Fig. [T]is due to a long zig-zag instability, which is not captured in our Landau ansatz. Similar 
remarks apply for instance to the rectangle branches, which are mostly unstable for Q over a 2 x 2 
domain. 



3.2 Ginzburg-Landau formalism and fronts 

Motivated by the acceptable approximation of the spots, stripes and beans via Q and (10) we 
proceed to use the Ginzburg-Landau (GL) reduction to predict the stationary fronts and localized 
patterns. If instead of Aj = Aj(t) we assume that Aj = Aj(t^x) are slowly varying functions also 
of X G M, then instead of the Landau system (10) we obtain the Ginzburg-Landau system 



dtAi = codlAi + h(Ai,A2,As), dtA, = JdlA, + f,{Ai,A2,A 



2,3, 



(13) 



where co(A) = -|a|^/x+((A:c, 0), A) > 0. See |Sch94[ lBvHS95l ISch99[ IMie02j for background on this 
formal procedure, and for so called at tr activity and approximation theorems which estimate the 
difference between a true solution of ^ and an approximation described by ([7| and ([13]), close to 
bifurcation. Thus, these theorems involve some small amplitude assumption for {u, v)^ related slow 
scales for t and x in Aj(t^x)^ and, for the present case of three resonant modes, a suitable scaling 
for the quadratic interactions, i.e., small C2. 
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Here we again want to use the GL system (13) at an 0(1) distance from Ac and find a necessary 
condition for stationary fronts between spots like A = (P+, P+, P+) and stripes like A = (T+, 0, 0) 



in (13). Thus we consider the stationary Ginzburg-Landau system 



codlAi+fi{Ai,A2,As) = 0, ^dlAj+fj{Ai,A2,A3) = 0, j = 2,3, 



(14) 



as a dynamical system in the spatial variable x. Now restricting to real amplitudes Aj^ the total 



energy of ([14|) is given by £;totai = ^kin + £^pot, where 



^ki] 



^0 //- . ^2 , 1/^ , x2 , 1/^ A\2 



2 Ad,A,y + -{d^A^r + -(d^As) 



and 



^POt - E (1^? + J^) + ^2^1^2A3 + |(^?^i + AlAl + Ai^: 

2=1 



^'^D 



are the kinetic and potential energy, respectively. Then ;^£^totai — 0, i.e., £^totai is conserved. Thus, 



and (T+,0,0) to 
Again we first focus on the hot bistable range 



a necessary condition for, e.g., a heteroclinic orbit Afront between (P^^P^^P^ 

exist in ( P i s £;pot(T+,0,0) = E^ot{P+.P+. P- 

and in Fig. 7(a) plot £'pot(7+, 0, 0) and £'pot(i^+, i^+, i^+). Their intersection defines the so called 

(hot) Maxwell point A^. Though we refrain from discussing the general energy landscape and 

dynamics of ( [m] ), it turns out that at least numerically the necessary condition A = A^ is also 

sufficient. Figure [TJ^b) shows a stationary front for (13) (with Neumann boundary conditions), 

which can either be obtained from time evolution of (13) with a suitable initial guess, or from 



solving the stationary boundary value problem. See also [MNT90] for some more analysis for front 
solutions of ( [m] ), including some implicit solution formulas. On the other hand, for A < A^ with 
Epot{T+i 0, 0) < £'totai(^+5 P+^ P+) we obtain a front for (13) travelling towards higher energy, i.e., 
(P+, P+, P+) invades (T+, 0, 0), and vice versa for A > A^. 
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Figure 7: (a), (c) Energies of spots and stripes near the hot resp. cold bistable regimes; their intersections 
define the respective Maxwell points A = 2.62, A = 3.11, and the zero of £^pot (P- , P- , P- ) defines the third 
Maxwell point A = 3.219. Hot (b) and cold (d) stationary Ginzburg-Landau fronts between stripes and 
spots at the respective Maxwell points, for comparison both over domain [— 127r//cc, 127r//cc], As = A2. 

Thus, from the GL approximation standing fronts t/het between different patterns are only 
predicted at precisely A = A^- In the physics literature, e.g., |Pom86j . the basic argument for the 
existence of t/het in an interval around A^ is that the patterns create an effective periodic potential 
which yields a pinning of fronts. 



It is not obvious whether the stationary GL system (14) has homoclinic solutions Ahom with, 
say, Ahom(^) -^ {P+^P+^P+) as X ^ ±oc, and that pass near (r+,0,0) near x = 0. However, 
(suitably shifted) fronts Afront(^o + •) and "backs" Aback(^i + •) with Aback(^) — ^front(~^) can be 
glued together to give approximate homoclinics with long plateaus but also with some dynamics in 
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time which can be expected to be exponentially slow in the separation distance between Afront ^ind 
^back? [CP90]. In fact we can generate almost stationary pulses numerically, but eventually the 



solution decays to a homogeneous rest state (T+,0,0) or {P^,P^,P^). Similarly, (14) may have 
periodic orbits which stay close to (P+, P+, P+) resp. (T^, 0, 0) over very long x-intervals. Clearly, 
for these "approximate homoclinics" a similar pinning argument as for the heteroclinics can be 
applied and predicts the existence of localized patterns as in Fig. ^h) near the Maxwell point A^. 
In [CKOOf [PMCKl 1 j it is worked out mathematically that the pinning and hence also the snaking 
are exponentially small effects and thus cannot be predicted at any order by Ginzburg-Landau type 
asymptotic expansions alone. However, combining a (high order) Ginzburg-Landau ansatz with 
beyond all order asymptotics, snaking in the model problems ([l]) and ^ can be described very 
accurately. See also, e.g., [K noOSf [ALB^IO and the references therein for alternative arguments 



explaining the snaking via so called heteroclinic tangles in the spatial dynamics formulation of 
again ^ and Q. 

The same pinning arguments apply to the cold bistable range, but the difference to the hot 
range is that cold Ginzburg-Landau fronts are flatter and wider since relatively to the vector fields 
/i, /2 the diffusion constant cq is much larger in the cold range than in the hot range. See Fig. [TJ^d). 
As a consequence, snaking might not exist between cold stripes and spots, or if it does exist, then 



it might be very difficult to find numerically over bounded domains, cf. Remark |2.5[ and similarly 
for fronts between the homogeneous solution and the cold hexagons. 

4 Discussion 

As for fronts between patterns and homogeneous solutions, the main ingredients for fronts t/het 
between two different periodic patterns is a bistable range, usually generated by a subcritical 
bifurcation (of mixed modes, in the case of two patterns), and the existence of a heteroclinic 
connection A^et between the corresponding fixed points in the associated Ginzburg-Landau system 
at some Maxwell point A = A^- The pinning effect in the full system then yields the existence of 
t/het in a parameter interval around A^, and the same effect yields branches of localized patterns 
Uhom in a patterned background even if the Ginzburg-Landau system only has "approximate" 
homoclinic solutions. 

These findings have already been suggested in [Pom86j , and been further worked out in ID 
problems and for patterns over homogeneous background in various papers, e.g., [B K07t [KnoOSt 
ICK09^ [ALB+10[ IDMCKll] . but the present paper appears to be the first numerical illustration for 
different patterns in 2D. Moreover, our numerics suggest that over bounded domains these branches 
may only snake if the bistability is sufficiently strong. These numerics are somewhat delicate due 
to the very many solution branches that exist in the Turing unstable range, in particular over large 
domains, which requires rather fine discretizations to avoid uncontrolled branch switching. 

Here we restricted to domains of still intermediate size, and this might be the reason why we 
were not able to generate snaking branches of fronts and localized patterns near the cold beans 
(Fig. |4]) or near Ac (Fig. [7|. For us it is an open roblem whether these branches start to snake 
over much larger domains. Close to the primary bifurcation we have a good Ginzburg-Landau 
approximation and possibly a beyond all orders expansion could be started |CK09^ IDMCKll] , but 
the snaking region might be very narrow and thus resolution of the snaking might require very 
large domains and very accurate numerics. 

The results are certainly not special for the model problem Q; we have used the same method 
to predict and numerically find various pinning branches in other 2D reaction diffusion systems 
|WF13j which are slightly more complicated algebraically due to Holling type nonlinearities, but 
for the same reason in fact nicer for the numerics. Interestingly, there all the pinning branches 
show snaking. 
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Remark 4.1. An ad hoc way to derive a Landau system like (10) from a system like Q is to only 
consider the solvability conditions (10) at the critical modes without first removing the A^ residual 
at second harmonics. This amounts to the ansatz w = J2i=i A^i^ + cc, or in other words setting 
005 01,02 = in ([9]). As this yields a residual of order A^ we call this the A^ ansatz and denote 
the new coefficients in (10) by C32 and C42. This is obviously simpler than ([t]), but not formally 
consistent. However, at order 0(1) distance from the bifurcation the A^ ansatz may give a better 
approximation of solutions than the residual of order A^ ansatz ([t]), since ^ represents only an 
asymptotic expansion, and not the first terms in some convergent series. 

It turns out that in some sense this is indeed the case for (Q, and this can be used to improve the 
prediction of the Maxwell point. While in the hot bistable range the A^ ansatz gives a much larger 
error for the stripes, its error for the hot spots is in fact smaller than the one in Fig. p[c). The idea 
is to use the coefficients cs{X) from the A^ ansatz for the stripes and the coefficients C32(A) from 
the A^ ansatz for the spots (|c4(A) — C42(A)| is very small anyway) in a "mixed" Ginzburg-Landau 
system that retains the variational structure. Thus, let S be T+ and H be P+ determined by using 
the A^- and A^-ansatz, respectively, and consider the system 

codlAi + ciAi + C2A2A3 + ( C33 ( ^l' " 1 + C32 ( ^i^ — ^ ) ) ^1 + c4Ai{Ai + Ai) = 0, 



4 

Co 



(Ai- 


-H 


U- 


-H 


fA2- 


-H 


\s- 


-H 


(A,- 


-H 



'^-dlA2 + Ci^2 + C2A1A3 + ( C33 ( ^-r I + C32 ( -^ I I ^1 + c^A2{A\ + Al) = 0, (15) 



""^-dlA-i + ci^3 + C2A1A2 + ( C33 ( ^1" " ) + C32 ( ^^ — ^ ) ) ^3 + 64^3(^1 + ^2) = 0. 



(Ar- 


-S 


\H- 


-S 


(A2- 


-s 


[h- 


-s 


fAs- 


-s 



In the hot bistable range we have S > H such that (14) is well-defined in this range. Moreover, 



(15) has again a conserved energy £'totai — ^kin + £'pot where now 



E,o^-j:'^A] + c,A,A2As + cssAt'^^-^ + cs2jZAti^^ + '^^^^^ 



Using this energy to calculate the Maxwell-point for (15) we find A^ = 2.67 which is more in the 
center of the snaking region than A^^ = 2.62. Thus, using the additional information that the spots 
are better approximated by the A^ ansatz than by the A^ ansatz we can obtain a better prediction 
for the Maxwell point. J 
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